Model fitting exercise

Step 1

Create a diagram named X and run the following commands to create a simple model configuration:

aadd "X" "POINT" "PO1" (100,110)
aadd "X" "POINT" "PO2" (140,110)
aadd "X" "POINT" "PO3" (180,100)
aadd "X" "POINT" "PO4" (180,120)
aadd "X" "PIPE" "PI1" (120,110)
aadd "X" "PIPE" "PI2" (160,100)
aadd "X" "PIPE" "PI3" (160,120)
amodi "PI1" "PI12_CONNECT_POINT_1" "PO1"
amodi "PI1" "PI12_CONNECT_POINT_2" "PO2"
amodi "PI2" "PI12_CONNECT_POINT_1" "PO2"
amodi "PI2" "PI12_CONNECT_POINT_2" "PO3"
amodi "PI3" "PI12_CONNECT_POINT_1" "PO2"
amodi "PI3" "PI12_CONNECT_POINT_2" "PO4"
aexclude "PO1"
aexclude "PO3"
aexclude "PO4"
amodi "PI1" "PI12_AREA" 0.06
amodi "PI2" "PI12_AREA" 0.03
amodi "PI3" "PI12_AREA" 0.03
amodi "ECCO" "MAXIMUM_TIME_STEP" 0.05
amodi "ECCO" "CURRENT_TIME_STEP" 0.05
amodi "SPEED" "SC_SPEED" 1000.0

Save the initial condition and set the simulation time to zero.

Step 2

Run the following commands to print how mass flow of the pipe PI1 behaves in a simple transient.

loadInitialCondition (syncRead $ \_ -> fromResource $ relativeResource currentModel "/Initial%20Condition")
runSequence mdo
      execute (setVar "PO1#PO11_PRESSURE" 1.1)
      fork (wait 1 >> stop)
      fork $ repeatForever mdo
          wait 0.05
          execute do
              massFlow = getVar "PI1#PI12_MIX_MASS_FLOW" :: Double
              print massFlow

Step 3

Modify the above commands so that they print the index of the time step and variables PI1#PI12_MIX_MASS_FLOW, PI2#PI12_MIX_MASS_FLOW and PI3#PI12_MIX_MASS_FLOW. The commands should produce something like:

0 3.4978646673069216 1.728445205173922 1.728445205173922
1 143.55517382976987 71.78198788751189 71.78198788751189
2 261.95050765748715 130.97625819564252 130.97625819564252
3 354.3858379588669 177.19356818048865 177.19356818048865
...

Use the following functions to maintain the counter (although it is possible to compute it also from the current simulation time).

Step 4

Modify the commands so that they compute the squared sum of errors from the "measurements" below:

measurements = [
    [0.05, 3.4977860762417605, 1.7283740169550754, 1.7284367287300313],
    [0.1, 134.88381983045093, 64.08353845423913, 70.80714606160574],
    [0.15, 231.9062548984293, 105.65792385470164, 126.24691709056889],
    [0.2, 296.4661748497294, 130.5064306373064, 165.95821168905044],
    [0.25, 336.98201115920403, 144.68635329210002, 192.29464808262938],
    [0.3, 361.56600347367186, 152.64673170010516, 208.9187479324867],
    [0.35, 376.0973084183755, 157.07188101607773, 219.0252113124207],
    [0.4, 384.690502383305, 159.57004084578205, 225.12040809039868],
    [0.45, 389.742925383871, 160.98719663407357, 228.75573737254442],
    [0.5, 392.7040735842458, 161.79511170873698, 230.90898615966157],
    [0.55, 394.4365833534546, 162.2577698637119, 232.17883313521224],
    [0.6, 395.36535200557034, 162.50160633728285, 232.86375430699712],
    [0.65, 395.93996807514526, 162.65081127152007, 233.28915819466746],
    [0.7, 396.2954901434386, 162.74229076331505, 233.55318938336018],
    [0.75, 396.51546801100426, 162.79847242004152, 233.71697886953123],
    [0.8, 396.65158867682203, 162.83302418223676, 233.81854447407886],
    [0.85, 396.7358254901475, 162.85429827014457, 233.881506261103],
    [0.9, 396.7879581657287, 162.8674097925786, 233.92052799957008],
    [0.95, 396.8202241096274, 162.87549715428165, 233.94470806962704],
    [1.0, 396.8401951361238, 162.88048891915273, 233.95968928259308]]

The columns of the data are, the simulation time PI1#PI12_MIX_MASS_FLOW, PI2#PI12_MIX_MASS_FLOW and PI3#PI12_MIX_MASS_FLOW.

Step 5

Now create a function objFun :: [Double] -> <Proc> Double that is called as

objFun [1.0,1.1,1.2]

It should run the same commands as in previous steps and return the squared sum of errors, but additionally set PI1#PI12_LOSS_COEFF, PI2#PI12_LOSS_COEFF and PI3#PI12_LOSS_COEFF to the given values.

It is good idea to also print the parameter values at the beginning of the function (and remove all other printing).

Step 6

Use the function

to fit the loss coefficients to the measurements. Try first with small number of function evaluations, because the current implementation does not allow interrupting the optimization. You may try rhobeg=1 and rhoend=1e-4 and initial guess [1,1,1].